module buty
  !
  ! Fortran Interface for butyA.s
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
  implicit none
  public mkhis_init
  public mkhisqh
!  public mkhisq Deprecated
  private
  character*(*),parameter::MEMFAIL="*** cannot allocate memory."
  character*(*),parameter::READFAIL="*** error during read."
  integer,parameter::MAXPATHLEN=1024
  integer,target::ndh=4           ! number of blocks in a plane
  integer nd                      ! nd*2
  integer,parameter::nt=64        ! number of crystals in transaxial
  integer,target::na=20           ! number of crystals in axial
  integer,target::mgap=4          ! number of crystals in a gap
  integer,target::nx=92           !na*ndh+mgap*(ndh-1) ! number of crystals in x
  integer,target::ny=nt           ! number of crystals in y
  integer nxy
  !
  integer*8,parameter::fbs=1024*1024*1024 ! maxmum byte to read at once
  integer*8,parameter::nemax=fbs/4
  integer*4,allocatable::buf(:)
  integer nbuf
  character*(MAXPATHLEN) fin,str
  character*(MAXPATHLEN),pointer::foutt,foutp,foutd,flist
  integer istat
  integer*8 n,ie,b1,b2,fs,ne
  integer ifile,nf,id,if1
  integer flg
  integer ix1,iy1,ix2,iy2
  logical lo_found
  integer*8 nb
  integer npair,ndelay
  logical::overwrite=.false.
  logical,allocatable::isfile(:)
  logical::ignoreerr=.false.
  logical::printcmdline=.false.  
  logical::useflist=.false.
  integer,pointer::indx(:)
  integer narg
  integer filetype
  logical doecho

contains
  
  !
  !
  subroutine echo_red(mess,code)
    character*(*),intent(in)::mess
    integer,intent(in),optional::code
    integer ec
    character*(*),parameter::swan_song="*** emergency escape from buty."
    write(9,*) trim(adjustl(mess))
!    call set_rstat(RSTAT_ERROR)
    ec=-1
    if(present(code)) ec=code
    call exit(ec)
  end subroutine echo_red

  !
  subroutine echo(m)
    character*(*),intent(in)::m
    if(.not.doecho) return
    write(*,10) trim(m)
10  format(x,"(bt)",x,a)
  end subroutine echo

  !
  subroutine mkhis_init(y,a,g,d,n,t,e,i,f)
    integer,intent(in)::y,a,g,d
    integer,intent(in)::n
    integer,intent(in)::t
    logical,intent(in)::e
    integer,intent(in),target,optional::i(n)
    character*(*),intent(in),optional::f
    filetype=t
    ndelay=0
    npair=0
    nb=0
    nf=0
    nbuf=0
    ny=y
    na=a
    mgap=g
    ndh=d
    nd=ndh*2
    nx=na*ndh+mgap*(ndh-1) ! number of crystals in x
    narg=n
    fin=""
    doecho=e
    nullify(indx)
    if(present(i)) indx=>i
    if(present(f)) then
       narg=1
       fin=f
    end if
    nxy=nx*ny
  end subroutine mkhis_init
  
  integer function get_next_file(ifile)
    use suty
    use kuty
    integer,intent(inout)::ifile
    character*512 mess
    ifile=ifile+1
    get_next_file=0 ! cycle
    if(associated(indx)) then
       if(ifile<=narg) then
          if(indx(ifile)==filetype) then
             call get_command_argument(ifile,fin)
          else
             return
          end if
       else
          get_next_file=-1
          return
       end if
!!$    else if(useflist) then
!!$       read(3,"(a)",iostat=istat) fin
!!$       if(istat/=0) then
!!$          get_next_file=-1
!!$          return
!!$       end if
!!$       if(len_trim(fin)==0) return
    else if(fin/="") then
    else
       get_next_file=-1 ! exit
       return
    end if
    fin=trim(adjustl(fin))     
    call echo("Opening file: "//trim(fin))
    open(unit=1,file=fin,status="old",access="stream",iostat=istat)
    if(istat/=0) then
       mess = "*** cannot open input file: code = "//trim(i2a(istat))
       if(ignoreerr) then
          call echo(mess)
          return
       else 
          call echo_red(mess)
       end if
    end if
    nf=nf+1
    fs=fsize(fin)
    call echo("size = "//trim(i2a(fs))//" byte")
    if(fs==0) then
       close(1)
       return
    end if
    get_next_file=1
  end function get_next_file
  
  subroutine print_stats(file,byte,ok,delay,&
       pair,flip,over,under,lost)
    use suty
    integer,intent(in)::file,byte,ok,delay
    integer,intent(in),optional::pair,flip,over,under,lost
    
    call echo("Number of files read:    "//trim(i2a(file)))
    call echo("Number of bytes read:    "//trim(i2a(byte)))
    if(present(pair))&
         call echo("Number of pair:          "//trim(i2a(pair)))
    call echo("Number of accept events: "//trim(i2a(ok)))
    call echo("Number of random events: "//trim(i2a(delay)))
    if(present(flip))&
         call echo("Number of fliped events: "//trim(i2a(flip)))
    if(present(over))&
         call echo("Number of overflow:      "//trim(i2a(over)))
    if(present(under))&
         call echo("Number of underflow:     "//trim(i2a(under)))
    if(present(lost))&
         call echo("Number of lost 4-byte:   "//trim(i2a(lost)))
  end subroutine print_stats
  
!!$  subroutine print_stats(file,byte,ok,delay)
!!$    use kuty
!!$    integer,intent(in)::file,byte,ok,delay
!!$    call echo("Number of files read:    "//trim(i2a(file)))
!!$    call echo("Number of bytes read:    "//trim(i2a(byte)))
!!$    call echo("Number of accept events: "//trim(i2a(ok)))
!!$    call echo("Number of random events: "//trim(i2a(delay)))
!!$    call echo("Delayed/Prompt (%):      "//trim(f2a(100*real(delay)/real(ok),1)))
!!$  end subroutine print_stats

  subroutine mkhis(plist,szlist,ppd)
    use kuty
    use suty
    integer,intent(in)::plist,szlist,ppd
    interface
       subroutine decd(pe,pd,tad1,tad2)
         integer*8,intent(in),value::pe
         integer*4,intent(out)::pd
         integer*4,intent(out)::tad1(3)
         integer*4,intent(out)::tad2(3)
       end subroutine decd
    end interface
    integer nok,noverflow,nflip,nunderflow,nlost
    integer*4,target::tad(3,2)
    integer*4 pd,tad_max(3,2),tad_min(3,2)
    integer*4,pointer::t1,a1,d1,t2,a2,d2
    integer,allocatable::off(:)
    integer pbuf
    integer at1,at2
    integer*4 dw
    pointer(pdw,dw)
    integer*2 c
    pointer(p,c)

    nok=0
    noverflow=0
    nflip=0
    nunderflow=0
    noverflow=0
    nlost=0
    t1=>tad(1,1); a1=>tad(2,1); d1=>tad(3,1)
    t2=>tad(1,2); a2=>tad(2,2); d2=>tad(3,2)              
    
    !
    allocate(off(ndh))
    do id=1,ndh
       off(id)=mgap*(id-1)
    end do
    
    ifile=0
    file: do
       if(plist==0) then
          select case(get_next_file(ifile))
          case(-1) 
             exit
          case(0)
             cycle
          end select
          
          n=fs/4
          if(n>nbuf) then
             if(nbuf>0) deallocate(buf)
             allocate(buf(n),stat=istat)
             if(istat/=0) call echo_red(MEMFAIL)
             nbuf=n
          end if
          
          b1=1
          read: do ie=1,n/fbs+1             
             b2=min(b1+fbs,n)
             read(1,iostat=istat) buf(b1:b2)    
             if(istat/=0) call echo_red(READFAIL)
             nb=nb+(b2-b1+1)
             b1=b2+1
          end do read
          close(1)
          pbuf=loc(buf)
       else             
          pbuf=plist
          n=szlist/4
       end if
       
       lo_found=.false.
       
       decode: do ie=1,n
          if(lo_found) then
             pdw=pbuf+(ie-1)*4
             flg=get_flg(dw)
             
             if(flg==0) then                   
                
                ! pair found
                lo_found=.false.
                npair=npair+1
                pdw=pdw-4
                call decd(pdw,pd,tad(1,1),tad(1,2))
                tad_max=max(tad,tad_max)
                tad_min=min(tad,tad_min)
                !
                if(t1>=nt.or.t2>=nt.or.a1>=na.or.a2>=na) then
                   call echo("*** WARNING: index overflow detected.")
                   call echo("ie,t1,a1,t2,a2="//trim(i2a(ie))//", "//trim(i2a(t1))//", "&
                        //trim(i2a(a1))//", "//trim(i2a(t2))//", "//trim(i2a(a2)))
                   noverflow=noverflow+1
                   cycle
                end if
                if(t1<0.or.t2<0.or.a1<0.or.a2<0) then
                   call echo("*** WARNING: index underflow detected.")
                   call echo("ie,t1,a1,t2,a2="//trim(i2a(ie))//", "//trim(i2a(t1))//", "&
                        //trim(i2a(a1))//", "//trim(i2a(t2))//", "//trim(i2a(a2)))
                   nunderflow=nunderflow+1
                   cycle
                end if

                tad=tad+1_4 ! <<<                
                if(d1<=ndh) then
                   if(d2<=ndh) cycle
                else
                   if(d2>=ndh+1) cycle
                   nflip=nflip+1
                   call echo("*** WARNING: index flip detected.")
                   call echo("d1,d2="//trim(i2a(d1))//", "//trim(i2a(d2)))
                   CYCLE ! <<< dont know what to do but for safty
!!$                 t1=>tad(1,2); a1=>tad(2,2); d1=>tad(3,2)
!!$                 t2=>tad(1,1); a2=>tad(2,1); d2=>tad(3,1)
                end if
                
                if(d1>ndh.or.d2>nd) then
                   call echo("*** WARNING: index overflow detected.")
                   call echo("d1,d2="//trim(i2a(d1))//", "//trim(i2a(d2)))
                   noverflow=noverflow+1
                end if
                if(d1<=0.or.d2<=ndh) then
                   call echo("*** WARNING: index underflow detected.")
                   call echo("d1,d2="//trim(i2a(d1))//", "//trim(i2a(d2)))
                   nunderflow=nunderflow+1
                end if
                !
                ix1=(d1-1)*na+(na-a1+1)+off(d1)
                iy1=t1
                ix2=(d2-(ndh+1))*na+a2+off(d2-ndh)
                iy2=t2
                !
                if(ix1>nx.or.ix2>nx.or.iy1>ny.or.iy2>ny) then
                   call echo("*** WARNING: index overflow detected.")
                   call echo("t1,a1,d1,t2,a2,d2="//trim(i2a(t1))//", "//trim(i2a(a1))//", "//trim(i2a(d1))//", "&
                        //trim(i2a(t2))//", "//trim(i2a(a2))//", "//trim(i2a(d2)))
                   call echo("ix1,iy1,ix2,iy2="//trim(i2a(ix1))//", "//trim(i2a(iy1))//", "//trim(i2a(ix2))//", "//trim(i2a(iy2)))
                   noverflow=noverflow+1
                   cycle
                end if
                if(ix1<=0.or.ix2<=0.or.iy1<=0.or.iy2<=0) then
                   call echo("*** WARNING: index overflow detected.")
                   call echo("t1,a1,d1,t2,a2,d2="//trim(i2a(t1))//", "//trim(i2a(a1))//", "//trim(i2a(d1))//", "&
                        //trim(i2a(t2))//", "//trim(i2a(a2))//", "//trim(i2a(d2)))
                   call echo("ix1,iy1,ix2,iy2="//trim(i2a(ix1))//", "//trim(i2a(iy1))//", "//trim(i2a(ix2))//", "//trim(i2a(iy2)))
                   nunderflow=nunderflow+1
                   cycle
                end if
                at1=ix1-1+nx*(iy1-1)
                at2=ix2-1+nx*(iy2-1)
                p=ppd+(at1+nxy*(at2+nxy*pd))*2
                c=c+1_2
                if(pd/=1) ndelay=ndelay+1
                nok=nok+1
             else
                nlost=nlost+1
                if(flg==1) then
                   lo_found=.true.
                else
                   lo_found=.false.
                end if
             end if
          else        
             if(get_flg(buf(ie))==1) then
                lo_found=.true.
             else
                nlost=nlost+1
             end if
          end if
       end do decode
    end do file
    
    call print_stats(nf,nb*4,nok,ndelay,&
         npair,nflip,noverflow,nunderflow,nlost)
    
  contains
    
    integer function get_flg(e)
      integer*4,intent(in)::e
      get_flg=ishft(e,-30)
    end function get_flg
    
  end subroutine mkhis
  
  !
  subroutine mkhisqh(plist,szlist,ppd)
    use suty
    interface
       subroutine initdecvqh(a,b,c,d,e)
         integer*8,intent(in),value::a,b,c,d,e
       end subroutine initdecvqh
    end interface
    interface
       integer*8 function decvqh(pe,ne,p,pnd)
         integer*8,intent(in),value::pe,ne,p,pnd
       end function decvqh
    end interface
    integer,intent(in)::plist,szlist,ppd
    integer pbuf,nread

    if(ndh/=4) then
       call mkhis(plist,szlist,ppd)
       return
       call echo_red("*** Number of blocks must be 4, use mkhis otherwise.")
    end if
    
    call initdecvqh(na,nt,ndh,mgap,nxy)
    if(istat/=0) call echo_red(MEMFAIL)
    pbuf=0
    ifile=0
    file: do
       if(plist==0) then
          select case(get_next_file(ifile))
          case(-1)
             exit
          case(0)
             cycle
          end select
          ne=fs/4
          n=min(ne,nemax)     
          if(nbuf<n) then
             if(nbuf/=0) deallocate(buf)
             allocate(buf(n),stat=istat)
             if(istat/=0) call echo_red(MEMFAIL)
             nbuf=n
          end if
       end if
       
       nread=0
       read: do        
          if(plist==0) then
             read(1,iostat=istat) buf(1:n)
             if(istat/=0) call echo_red(READFAIL)
             nread=nread+n
             nb=nb+n
             pbuf=loc(buf)
          else
             pbuf=plist
             n=szlist/4
          end if
          
          npair=npair+decvqh(pbuf,n,ppd,loc(ndelay))
          
          if(plist/=0) exit file
          if(nread==ne) exit          
          n=min(ne-nread,nemax)
          
       end do read
       
       close(1) ! <<<
       
    end do file
    
    call print_stats(nf,nb*4,npair,ndelay)
    
  end subroutine mkhisqh

!
!!$ Deprecated
!
!!$  subroutine mkhisq(plist,szlist,ppd)
!!$    use suty
!!$    interface
!!$       subroutine initdecvq(a,b,c,d)
!!$         integer*8,intent(in),value::a,b,c,d
!!$       end subroutine initdecvq
!!$    end interface
!!$    interface
!!$       subroutine decvq(pe,ne,pb,nb,pneo,pnbo)
!!$         integer*8,intent(in),value::pe,ne,pb,nb,pneo,pnbo
!!$       end subroutine decvq
!!$    end interface
!!$    integer,intent(in)::plist,szlist,ppd
!!$    integer*8,parameter::nlmdbufmax=737280
!!$    integer nlmdbuf,m
!!$    integer*4,allocatable,target::lmdbuf(:,:)
!!$    integer*8 nlmdbuf_wrote, nbuf_read,nread
!!$    integer pbuf
!!$    integer*4,pointer::at1,at2,pd
!!$    integer*2 c
!!$    pointer(p,c) 
!!$    if(ndh/=4) then
!!$       call mkhis(plist,szlist,ppd)
!!$       return
!!$!       stop "*** Number of blocks must be 4, use mkhis otherwise."
!!$    end if
!!$    call initdecvq(na,nt,ndh,mgap)
!!$    if(istat/=0) stop MEMFAIL   
!!$    pbuf=0
!!$    nlmdbuf=0    
!!$    ifile=0
!!$    file: do
!!$       if(plist==0) then
!!$          select case(get_next_file(ifile))
!!$          case(-1)
!!$             exit
!!$          case(0)
!!$             cycle
!!$          end select
!!$          
!!$          ne=fs/4
!!$          n=min(ne,nemax)     
!!$          if(nbuf<n) then
!!$             if(nbuf/=0) deallocate(buf)
!!$             allocate(buf(n),stat=istat)
!!$             if(istat/=0) stop MEMFAIL
!!$             nbuf=n
!!$          end if
!!$          
!!$          m=min(fs/8,nlmdbufmax)
!!$          if(m>nlmdbuf) then
!!$             if(nlmdbuf/=0) deallocate(lmdbuf)
!!$             allocate(lmdbuf(4,m),stat=istat)
!!$             if(istat/=0) stop MEMFAIL
!!$             nlmdbuf=m
!!$          end if
!!$       end if
!!$       
!!$       nread=0
!!$       read: do        
!!$          if(plist==0) then
!!$             read(1,iostat=istat) buf(1:n)
!!$             if(istat/=0) stop READFAIL
!!$             nread=nread+n
!!$             nb=nb+n
!!$             pbuf=loc(buf)
!!$          else
!!$             pbuf=plist
!!$             n=szlist/4
!!$          end if
!!$          
!!$          dec: do
!!$             
!!$             call decvq(pbuf,n,loc(lmdbuf),nlmdbufmax,&
!!$                  loc(nbuf_read),loc(nlmdbuf_wrote))
!!$             npair=npair+nlmdbuf_wrote
!!$             
!!$             sort: do ie=1,nlmdbuf_wrote
!!$                at1=>lmdbuf(1,ie); at2=>lmdbuf(2,ie)
!!$                pd=>lmdbuf(3,ie)
!!$                if(pd==0) ndelay=ndelay+1
!!$                p=ppd+(at1+nxy*(at2+nxy*pd))*2
!!$                c=c+1
!!$             end do sort
!!$             
!!$             if(nbuf_read==n) exit
!!$             pbuf=pbuf+4*nbuf_read
!!$             n=n-nbuf_read
!!$             
!!$          end do dec
!!$
!!$          if(plist/=0) exit file
!!$          if(nread==ne) exit          
!!$          n=min(ne-nread,nemax)
!!$          
!!$       end do read
!!$       
!!$       close(1) ! <<<
!!$       
!!$    end do file
!!$    
!!$    call print_stats(nf,nb*4,npair,ndelay)
!!$
!!$  end subroutine mkhisq
  
end module buty
